rm(list = ls())
# preparing data ------------------------------------------------------------
library(tidyverse)
library(ggthemes)
library(grid)
library(gridExtra)
library(ggpubr)
library(lubridate)
library(estimatr)
library(texreg)
uk = readRDS("ukraine_aid.rds")


# Generating levels for demographic factors
uk = uk %>% 
  mutate(age2 = cut(age, quantile(age, c(0, .33, .67, 1), na.rm = T), include.lowest = T, labels = c("18-30", "31-40", "41-80")),
         educ2 = cut(educ, c(0, 6, 7), include.lowest = T, labels = c("Less than University", "University or Higher")),
         income2 = cut(income, breaks = quantile(income, c(0, 0.33, 0.67, 1), na.rm = T), include.lowest = T, labels = c("< 5000", "5000-8000", "> 8000")),
         ethnicity2 = recode(ethnicity, "I feel more Russian than Ukrainian" = "Russian", "I feel only Russian" = "Russian",
                             "I feel more Ukrainian than Russian" = "Ukrainian", "I feel only Ukrainian" = "Ukrainian", 
                             "I feel equally Ukrainian and Russian" = "Both", "I feel neither Russian nor Ukrainian" = "Neither"),
         language2 = recode(language, "I speak both, but mostly Russian" = "Russian", "I speak Russian" = "Russian",
                            "I speak both, but mostly Ukrainian" = "Ukrainian", "I speak Ukrainian" = "Ukrainian", 
                            "I speak both equally" = "Both", "I speak surzhik" = "Surzhik"),
         job2 = recode(job, "Military servant" = "Other", "Entrepreneur, farmer" = "Other", "Pension (because of age or disability)" = "Other",
                       "Student" = "Other", "Worker, farmworker" = "Other", "Householder" = "Other", "Servant (without higher education)" = "Other",
                       " Other" = "Other", "Self employed businesswomen/men" = "Professional", 
                       "Professional (with higher education)" = "Professional"),
         job2 = factor(job2, levels = c("Other", "Professional", "Unemployed")),
         sex2 = factor(sex, levels = c("Male", "Female")),
         eu_view2 = recode(eu_view, "0" = "Unfavorable", "1" = "Unfavorable", "2" = "Neutral", 
                           "3" = "Favorable", "4" = "Favorable"),
         eu_view2 = factor(eu_view2, levels = c("Unfavorable", "Neutral", "Favorable")),
         rus_view2 = recode(rus_view, "0" = "Unfavorable", "1" = "Unfavorable", "2" = "Neutral", 
                            "3" = "Favorable", "4" = "Favorable"),
         rus_view2 = factor(rus_view2, levels = c("Unfavorable", "Neutral", "Favorable")),
         eu_aid2 = recode(eu_aid, "0" = "Unfavorable", "1" = "Unfavorable", "2" = "Neutral", 
                          "3" = "Favorable", "4" = "Favorable"),
         eu_aid2 = factor(eu_aid2, levels = c("Unfavorable", "Neutral", "Favorable")),
         rus_aid2 = recode(rus_aid, "0" = "Unfavorable", "1" = "Unfavorable", "2" = "Neutral", 
                           "3" = "Favorable", "4" = "Favorable"),
         rus_aid2 = factor(rus_aid2, levels = c("Unfavorable", "Neutral", "Favorable"))
  )


# Descriptive results -----------------------------------------------------
p1 = uk %>% 
  filter(complete.cases(eu_view2)) %>% 
  ggplot() + 
  geom_bar(aes(x = eu_view2, y = ..count../sum(..count..)), fill = "gray60") + 
  theme_few() + 
  labs(x = "", title = "Views toward the EU", y = "Percent of Respondents") + 
  scale_y_continuous(labels = scales::percent) + 
  coord_cartesian(ylim = c(0, 0.4))

p2 = uk %>% 
  filter(complete.cases(rus_view2)) %>% 
  ggplot() + 
  geom_bar(aes(x = rus_view2, y = ..count../sum(..count..)), fill = "gray60") + 
  theme_few() + 
  labs(x = "", title = "Views toward Russia", y = "Percent of Respondents") +
  scale_y_continuous(labels = scales::percent) + 
  coord_cartesian(ylim = c(0, 0.4))

ggarrange(p1, p2, ncol = 2) %>% ggexport(filename = "figure1.pdf", width = 7, height = 5)

p1 = uk %>% 
  filter(complete.cases(eu_aid2)) %>% 
  ggplot() + 
  geom_bar(aes(x = eu_aid2, y = ..count../sum(..count..)), fill = "gray60") + 
  theme_few() + 
  labs(x = "", title = "Views toward EU Aid", y = "Percent of Respondents") + 
  scale_y_continuous(labels = scales::percent) + 
  coord_cartesian(ylim = c(0, 0.6))

p2 = uk %>% 
  filter(complete.cases(rus_aid2)) %>% 
  ggplot() + 
  geom_bar(aes(x = rus_aid2, y = ..count../sum(..count..)), fill = "gray60") + 
  theme_few() + 
  labs(x = "", title = "Views toward Russian Aid", y = "Percent of Respondents") + 
  scale_y_continuous(labels = scales::percent) + 
  coord_cartesian(ylim = c(0, 0.6))

ggarrange(p1, p2, ncol = 2) %>% ggexport(filename = "figure2.pdf", width = 7, height = 5)

# Framing ------------------------------------------------------------
#mean favorability of aid
uk = uk %>% 
  rowwise() %>% 
  mutate(aid_mean = mean(c(eu_aid, rus_aid)))

mod1 = lm_robust(eu_aid ~ t_aid, data = uk)

mod2 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk)

mod3 = lm_robust(rus_aid ~ t_aid, data = uk)

mod4 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk)

#capture and save output from texreg
texreg(list(mod1, mod2, mod3, mod4), include.ci = F, 
       omit.coef = "age|educ|income|ethnicity|language|job|sex|region|pcr|fam",
       custom.coef.names = c("Constant", "Humanitarian Frame", "Political Influence Frame"),
       label = "tab:framing",
       caption = "Treatment effects on attitudes towards aid",
       include.rmse = F, reorder.coef = c(3, 2, 1), stars = c(0.01, 0.05, 0.1))

# Favorable ---------------------------------------------------------------
mod1 = lm_robust(eu_aid ~ t_aid, data = uk[uk$eu_view>2,])
mod2 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view>2,])
mod3 = lm_robust(rus_aid ~ t_aid, data = uk[uk$rus_view>2,])
mod4 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$rus_view>2,])

texreg(list(mod1, mod2, mod3, mod4), include.ci = F, 
                            omit.coef = "age|educ|income|ethnicity|language|job|sex|region|pcr|fam",
                            custom.coef.names = c("Constant", "Humanitarian Frame", "Political Influence Frame"),
                            label = "tab:favorable", float.pos = "H",
                            caption = "Treatment effects on attitudes towards aid among those favorable of foreign actors",
                            include.rmse = F, reorder.coef = c(3, 2, 1), stars = c(0.01, 0.05, 0.1))


# Neutral -----------------------------------------------------------------
mod1 = lm_robust(eu_aid ~ t_aid, data = uk[uk$eu_view==2,])
mod2 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view==2,])
mod3 = lm_robust(rus_aid ~ t_aid, data = uk[uk$rus_view==2,])
mod4 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$rus_view==2,])

texreg(list(mod1, mod2, mod3, mod4), include.ci = F, 
                            omit.coef = "age|educ|income|ethnicity|language|job|sex|region|pcr|fam",
                            custom.coef.names = c("Constant", "Humanitarian Frame", "Political Influence Frame"),
                            label = "tab:neutral", float.pos = "H",
                            caption = "Treatment effects on attitudes towards aid among those neutral toward foreign actors",
                            include.rmse = F, reorder.coef = c(3, 2, 1), stars = c(0.01, 0.05, 0.1))

# Unfavorable -----------------------------------------------------------
mod1 = lm_robust(eu_aid ~ t_aid, data = uk[uk$eu_view<2,])
mod2 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view<2,])
mod3 = lm_robust(rus_aid ~ t_aid, data = uk[uk$rus_view<2,])
mod4 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$rus_view<2,])

texreg(list(mod1, mod2, mod3, mod4), include.ci = F, 
                            omit.coef = "age|educ|income|ethnicity|language|job|sex|region|pcr|fam",
                            custom.coef.names = c("Constant", "Humanitarian Frame", "Political Influence Frame"),
                            label = "tab:unfavorable", float.pos = "H",
                            caption = "Treatment effects on attitudes towards aid among those unfavorable of foreign actors",
                            include.rmse = F, reorder.coef = c(3, 2, 1), stars = c(0.01, 0.05, 0.1))


# Discussion figure -------------------------------------------------------
#Favorable
mod1 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view>2,])
mod2 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$rus_view>2,])

#Neutral
mod3 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view==2,])
mod4 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$rus_view==2,])


#Unfavorable
mod5 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view<2,])
mod6 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$rus_view<2,])

df_plot = tidy(mod1) %>% mutate(view = "Favorable") %>% bind_rows(tidy(mod2) %>% mutate(view = "Favorable")) %>% 
  bind_rows(tidy(mod3) %>% mutate(view = "Neutral")) %>% bind_rows(tidy(mod4) %>% mutate(view = "Neutral")) %>% 
  bind_rows(tidy(mod5) %>% mutate(view = "Unfavorable")) %>% bind_rows(tidy(mod6) %>% mutate(view = "Unfavorable")) %>% 
  filter(str_detect(term, "^t_")) %>% 
  mutate(view = factor(view, levels = c("Unfavorable", "Neutral", "Favorable")))

df_plot %>% 
  mutate(term = recode(term, "t_aidhumanitarian" = "Humanitarian Frame", "t_aidpolitical" = "Political Influence Frame"),
         outcome = recode(outcome, "eu_aid" = "EU Aid", "rus_aid" = "Russian Aid")) %>% 
  ggplot(aes(x = view, y = estimate, ymin = conf.low, ymax = conf.high, position = term, colour = term)) + 
  geom_point(position = position_dodge(width = .5)) + 
  geom_pointrange(position = position_dodge(width = .5)) + 
  geom_hline(yintercept = 0, linetype = 2) + 
  coord_flip() + 
  facet_wrap(~outcome) + 
  theme_few() + 
  theme(legend.position = "bottom", legend.title = element_blank()) + 
  labs(x = "", y = "") + 
  scale_color_grey(start = .2, end = .7)
ggsave("figure3.pdf", width = 5, height = 5)

# East out --------------------------------------------------------------------
mod1 = lm_robust(eu_aid ~ t_aid, data = uk[uk$region == "east",])
mod2 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + pcr_mod, data = uk[uk$region == "east",])
mod3 = lm_robust(rus_aid ~ t_aid, data = uk[uk$region == "east",])
mod4 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + pcr_mod, data = uk[uk$region == "east",])

texreg(list(mod1, mod2, mod3, mod4), include.ci = F, 
                            omit.coef = "age|educ|income|ethnicity|language|job|sex|region|pcr|fam",
                            custom.coef.names = c("Constant", "Humanitarian Frame", "Political Influence Frame"),
                            label = "tab:east", float.pos = "H",
                            caption = "Treatment effects on attitudes towards aid among those residing in the East",
                            include.rmse = F, reorder.coef = c(3, 2, 1), stars = c(0.01, 0.05, 0.1))


# sep --------------------------------------------------------------------
mod1 = lm_robust(eu_aid ~ t_aid, data = uk[uk$region == "sep",])
mod2 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + pcr_mod, data = uk[uk$region == "sep",])
mod3 = lm_robust(rus_aid ~ t_aid, data = uk[uk$region == "sep",])
mod4 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + pcr_mod, data = uk[uk$region == "sep",])

texreg(list(mod1, mod2, mod3, mod4), include.ci = F, 
                            omit.coef = "age|educ|income|ethnicity|language|job|sex|region|pcr|fam",
                            custom.coef.names = c("Constant", "Humanitarian Frame", "Political Influence Frame"),
                            label = "tab:sep", float.pos = "H",
                            caption = "Treatment effects on attitudes towards aid among those residing in the Separatist Republics",
                            include.rmse = F, reorder.coef = c(3, 2, 1), stars = c(0.01, 0.05, 0.1))



# Subset to EU favorable and Russia unfavorable --------------------------------------------
mod1 = lm_robust(eu_aid ~ t_aid, data = uk[uk$eu_view>2 & uk$rus_view < 2,])
mod2 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view>2 & uk$rus_view < 2,])
mod3 = lm_robust(rus_aid ~ t_aid, data = uk[uk$eu_view>2 & uk$rus_view < 2,])
mod4 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view>2 & uk$rus_view < 2,])

texreg(list(mod1, mod2, mod3, mod4), include.ci = F, 
                            omit.coef = "age|educ|income|ethnicity|language|job|sex|region|pcr|fam",
                            custom.coef.names = c("Constant", "Humanitarian Frame", "Political Influence Frame"),
                            label = "tab:eu_fav_ru_unfav", float.pos = "H",
                            caption = "Treatment effects on attitudes towards aid among those favorable of the EU and unfavorable of Russia",
                            include.rmse = F, reorder.coef = c(3, 2, 1), stars = c(0.01, 0.05, 0.1))


# Subset to EU unfavorable and Russia favorable --------------------------------------------
mod1 = lm_robust(eu_aid ~ t_aid, data = uk[uk$eu_view<2 & uk$rus_view > 2,])
mod2 = lm_robust(eu_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view<2 & uk$rus_view > 2,])
mod3 = lm_robust(rus_aid ~ t_aid, data = uk[uk$eu_view<2 & uk$rus_view > 2,])
mod4 = lm_robust(rus_aid ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk[uk$eu_view<2 & uk$rus_view > 2,])

texreg(list(mod1, mod2, mod3, mod4), include.ci = F, 
                            omit.coef = "age|educ|income|ethnicity|language|job|sex|region|pcr|fam",
                            custom.coef.names = c("Constant", "Humanitarian Frame", "Political Influence Frame"),
                            label = "tab:eu_unfav_ru_fav", float.pos = "H",
                            caption = "Treatment effects on attitudes towards aid among those favorable of Russia and unfavorable of the EU",
                            include.rmse = F, reorder.coef = c(3, 2, 1), stars = c(0.01, 0.05, 0.1))


# Ordinal Logistic regressions-------------------------------------------------------
library(MASS)
mod1 = polr(eu_aid2 ~ t_aid, data = uk, Hess = TRUE)
mod2 = polr(eu_aid2 ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk, Hess = TRUE)
mod3 = polr(rus_aid2 ~ t_aid, data = uk, Hess = TRUE)
mod4 = polr(rus_aid2 ~ t_aid + age + educ + income + ethnicity2 + language2 + job2 + sex2 + region + pcr_mod, data = uk, Hess = TRUE)

stargazer::stargazer(mod1, mod1, mod2, mod2, mod3, mod3, mod4, mod4,
          coef = list(coef(mod1), exp(coef(mod1)),
                      coef(mod2), exp(coef(mod2)),
                      coef(mod3), exp(coef(mod3)),
                      coef(mod4), exp(coef(mod4))), 
          se = list(sqrt(diag(vcov(mod1)))[1:2], c(10, 10),
                    sqrt(diag(vcov(mod2)))[1:2], c(10, 10),
                    sqrt(diag(vcov(mod3)))[1:2], c(10, 10),
                    sqrt(diag(vcov(mod4)))[1:2], c(10, 10)), dep.var.labels.include = F,
          column.labels = c("EU Aid", "EU (OR)", "EU", "EU (OR)",
                          "Russia", "Russia (OR)", "Russia", "Russia (OR)"), 
          omit = "age|educ|income|ethnicity|language|job|sex|region|pcr|fam", 
          covariate.labels = c("Humanitarian Frame", "Political Influence"),
          add.lines = list(c("Controls", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes")), 
          label = "tab:ordinal", table.placement = "H", column.sep.width = "0pt",
          title = "Treatment effects on attitudes towards aid using ordinal logistic regression", font.size = "footnotesize")


# Balance tables ----------------------------------------------------------
library(arsenal)
uk$gender = uk$sex
levels(uk$gender) = rev(levels(uk$gender))

uk$income2 = cut(uk$income, 
                 breaks = quantile(uk$income, probs = c(0, .33, .67, 1), na.rm = T), 
                 labels = c("Up ot 4000 hryvnia (33rd percentile)", "Up ot 7000 hryvnia (67th percentile)", "Over 7000 hryvnia"))

uk$ethnicity2 = str_replace(uk$ethnicity, "I feel ", "")
uk$language2 = str_replace(uk$language, "I speak", "Speak")
uk$religion2 = recode(uk$religion, "Greek Catholic Church" = "Other", "Muslim (Islam)" = "Other",
                      "Protestant" = "Other", "Roman Catholic" = "Other", 
                      "Ukrainian Autocephalous Orthodox Church" = "Other",
                      "Religious but do not belong to a certain religion or church" = "Religious but no church")

uk$region2 = recode(uk$region, "east" = "East", "sep" = "Separatist")

uk$job2 = recode(uk$job, "Entrepreneur, farmer" = "Other", "Military servant" = "Other",
                 "Pension (because of age or disability)" = "Pension", 
                 "Professional (with higher education)" = "Professional", 
                 "Self employed businesswomen/men" = "Self-employed",
                 "Servant (without higher education)" = "Servant (no higher ed.)")

uk$t_aid = recode(uk$t_aid, "control" = "Control", "humanitarian" = "Humanitarian", "political" = "Political Influence")
uk$university = ifelse(uk$educ == 7, 1,
                       ifelse(uk$educ < 7, 0, NA))

vars1 = c("age", "income2", "ethnicity2", "language2", "gender")
vars2 = c("religion2", "region2", "job2", "university")

formula1 = paste0("t_aid", "~", paste0(vars1, collapse = "+")) %>% as.formula()
formula2 = paste0("t_aid", "~", paste0(vars2, collapse = "+")) %>% as.formula()

tab1 = tableby(formula1, data = uk, total = F, digits = 2L, 
               control = tableby.control(numeric.stats = c("Nmiss2", "meansd", "range"),cat.stats = c("Nmiss2", "countpct"))) %>% 
  summary(labelTranslations = list(age = "Age (Years)", income2 = "Income", ethnicity2 = "Ethnicity",
                                   language2 = "Language", gender = "Gender"),
          pfootnote = F, text = T)

tab2 = tableby(formula2, data = uk, total = F, digits = 2L, 
               control = tableby.control(numeric.stats = c("Nmiss2", "meansd", "range"),cat.stats = c("Nmiss2", "countpct"))) %>% 
  summary(labelTranslations = list(religion2 = "Religion", region2 = "Region", job2 = "Job", 
                                   university = "University Educated"),
          pfootnote = F, text = T)

tab1 %>% as.data.frame() %>% xtable::xtable(caption = "Balance Table") %>% 
  xtable::print.xtable(include.rownames = F, floating = T, table.placement = "H", size = "footnotesize")

tab2 %>% as.data.frame() %>% xtable::xtable(caption = "Balance Table") %>% 
  xtable::print.xtable(include.rownames = F, floating = T, table.placement = "H", size = "scriptsize")
